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ABSTRACT 

The measurement of cosmological parameters is investigated in a representation of the 
least-action method that uses a redshift-space dataset to simultaneously constrain the 
real-space fields 5,v. This method is robust in recovering the entire evolution of the 
matter density contrast and peculiar velocities of galaxies in real space from current 
galaxy redshift surveys. The main strength of the method is that it permits us to break 
the degeneracy of the parameters b and Q m (customarily measured in the ratio j3 = 
0"; 6 /6 from redshift-space distortions), and these are evaluated in the current context 
separately. The procedure provides a simple numerical means to extract as much 
information as possible from a given sample, in the simplest linear bias model, before 
resorting to cosmic complementarity to resolve the degeneracy in the measurement 
of O m . The same premise applies to more sophisticated choices of bias models. We 
construct a likelihood parameter X(b, f2 m ) to evaluate the relative likelihood of different 
values of b and f2 m . The method is applied to the IRAS 1.2 Jy redshift survey with 
a low-resolution Gaussian smoothing length of 1200 km s _1 within a spherical region 
Xjuax ~ 15,000 km s _1 and the reconstructed velocity field is then compared with 
POTENT-reconstructed velocities from the Mark III radial-velocity dataset within a 
radius ~ 5000 km s~ , which have been suitably prepared to account for Malmquist 
bias and other systematic errors. The analysis yields a likelihood for the parameters 
that is overall consistent with Q m ks 0.3 and b ks 1.1, thus lending support to a 
non-vanishing cosmological constant w 0.7 in a flat universe. 

Key words: large-scale structure of universe - cosmology - galaxies: distances and 
redshifts 



1 INTRODUCTION 

Galaxy redshift surveys are undoubtedly extremely valuable 
tools to investigate the evolution of the universe at large 
scales. The cosmologist's prerogative is to determine the evo- 
lution of the matter density contrast 8 and peculiar velocity 
v that yields such cosmic structure, customarily assuming 
that it formed solely by gravity, and the cosmological pa- 
rameters that determine their dynamics. In the standard 
paradigm of a FRW expanding universe, the interplay of 
both fields is governed by the density parameter fio and the 
Hubble parameter Ho- On the other hand, a relationship 
between the fields S,v and the survey data is established by 
adopting a bias model that purports the correlation between 
the z-space galaxy number-count and the underlying matter 
field. Devoid of such a relationship, the edifice of measuring 
cosmological parameters from galaxy redshift surveys has no 
foundation whatsoever. A standard working hypothesis, that 
I shall accept throughout this paper, is that of linear bias, 



i.e. b 2 = P (fc) g a i s / Pik) ma ttor (more elaborate bias models 
are propounded in e.g. Dekel & Lahav 1999). For simplicity 
we shall also leave out the scale-dependence of b. Therefore, 
three relevant parameters that are interesting to pin down 
from redshift-space samples are in this context Q m ,Ho and 
b. In this paper I shall be chiefly concerned with Q m and b 
(Ho will be scaled out with distance). 

Tracing back in time the matter fields takes us to an ini- 
tial epoch of fluctuations of very small amplitude 5 & 10 -4 , 
seeded by a period of inflationary expansion. At that point 
the information derived from the galaxy surveys connects 
with early-universe data such as the spectrum of fluctua- 
tions on the CMB. If the matter fields could realistically be 
traced back to such a primordial stage by integrating the 
equations of gravitational instability, then the statistics of 
the 5 field would be a potentially key discriminant to rule 
out cosmological models. For instance, non-gaussianity in 
the initial 5 field rules out most inflationary models, and 
only those leading to a non-Gaussian primordial spectrum 
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Figure 1. Qualitative distribution of errors in the reconstruction of the matter fields from redshift surveys. The error eo m the current 
redshift sample increases monotonically in (a) as the pcrturbativc solutions propagate eo to increasing amplitudes when integrated back 
in time; in the boundary condition problem of the LAP method, shown in (b), errors fluctuate between the fixed end-points. 



remain acceptable (such models are suggested in e.g. Linde, 
Sasaki & Tanaka 1999). 

Kaiser (1987) proposed measuring cosmological param- 
eters from redshift-space distortions by virtue of the fact 
that overdense regions appear to be flatter along the line- 
of-sight in redshift space. This distortion, quantified by the 
parameter (3 = fij; 6 /b, permits us to solve the equations 
for S,v, at least perturbatively (see e.g. Dekel 1994; Coles 
& Sahni 1995), and measurements of /3 have been investi- 
gated in much detail in the literature (Strauss & Willick 
1995; Dekel 1994,1999a; Dekel, Burstein & White 1997). 
Also, in view of the fact that the bias parameter is almost 
certainly dependent on the selected sample, estimates have 
been computed for Piras given bi for IRA S galaxies (Dekel 
et al. 1993; Fisher et al. 1995a; Willick et al. 1997a,b; Sigad 
et al. 1998; more recently from the PSCz sample, Canavezes 
et al. 1998; Tadros et al. 1999; Saunders et al. 2000) and from 
the Optical Redshift Survey (ORS) (Hudson et al. 1995; 
Santiago et al. 1995; Baker et al. 1998). The Mark III pe- 
culiar velocity survey similarly yields estimates of /3 from 
redshift distortions (Willick et al. 1995, 1996, 1997a,b; Dekel, 
Burstein & White 1997; Sigad et al. 1998). It is only beyond 
the linear approximation (i.e. S oc V-v) and, indeed, beyond 
the assumption of linear bias, that one can break down the 
degeneracy between f2 m and 6 and estimate these param- 
eters separately, rather than via /3 (Fry 1994; Bernardeau 
et al. 1995). Verde et al. (1998) achieved this by proposing 
the bispectrum as a measure of cosmological parameters, 
in a model of non-linear bias. In this paper we also pursue 
breaking the degeneracy of fi m and b from the redshift-space 
data and show that by using the least-action framework it 
is indeed possible to do so within the linear bias model. 

The least-action principle (LAP) was first used in the 
Local Group by Peebles (1989,1990). The trajectories of 
nearby galaxies were computed subject to two boundary 
conditions: vanishing initial velocities and fixed present po- 
sitions. This simple scenario of self-gravitating point-like 



masses with two boundary conditions produced an estimate 
of £2 m by fitting to the observations the predicted peculiar 
velocities of nearby galaxies. The LAP method has also been 
used as a test of f2 m = 1 CDM models (Branchini & Carl- 
berg 1994), as well as to integrate the orbits of a significant 
number of galaxies from partial coverage redshift samples 
(e.g. Shaya, Peebles & Tully 1995). An equivalent repre- 
sentation of the LAP method in terms of continuous fields, 
i.e. the density contrast and velocity fields was proposed 
by Giavalisco et al. (1993), and employed in Susperregi & 
Binney (1994) (hereafter SB94) and Susperregi (1995) in the 
reconstruction of fl m = 1 simple models, such as exact solu- 
tions and Gaussian random fields. More recently, Schmoldt 
& Saha (1998) proposed a variant of the customary LAP 
formulation by rewriting the equations motion in redshift 
space. 

The key difference between the variational and pertur- 
bative approaches lies on how the errors are spread over 
the time-reversed evolution. This is qualitatively sketched 
in Fig. 1. A nth-order solution differs, in the time-reversed 
direction, from the true solution by a monotonically grow- 
ing parameter e which sets out from a small value e(io) (at 
any rate eo is at least the sum of the systematic and random 
errors of the dataset) and the conservation of kinematical 
quantities is preserved up to 0(e n ). This is adequate within 
a time span t c <C t Ss to where e(£ c ) ~ 1, and t c marks 
a transition into the loss of convergence. The distribution 
of errors in the LAP method on the other hand, is by con- 
struction evenly distributed along the trajectory; the initial 
and final boundary conditions are fixed, though not without 
systematic and numerical errors, and the parameter e fluctu- 
ates along the trajectory between both end-points (Fig. lb). 
Hence the solution is well-behaved whether the errors remain 
within the bound e Ss 1 or not. In that respect there is an 
advantage with respect to perturbative solutions; the down- 
side of it is of course that within the span of time where 
perturbative solutions are valid, LAP errors may fluctuate 
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with larger amplitude than the perturbative equivalent. The 
LAP method, in a nutshell, thus consists in finding Ansatze 
for the matter fields that optimize the distribution of e along 
the phase-space trajectory, and hence minimize the overall 
departure with respect to the exact solution. The following 
two difficulties may arise: 

• A Finding "dynamically plausible" solutions. If the 
matter field is sparsely sampled or the errors in the dataset 
are substantial, then the boundary condition given by the 
survey, taken at face value, may not correspond to the out- 
come of gravitational evolution from the initial fluctuations 
(typically S ~ or vanishing peculiar velocities). The LAP 
method will in this case find a dynamically plausible fit be- 
tween the end-points, which will be as faithful a representa- 
tion of the true evolution as is the quality of the dataset. 

• B Formation of multistreams in over-dense regions. 
Multistreams are characterised by galaxies at the same red- 
shift which are located at different positions along the line 
of sight and have different infalling velocities. The degener- 
acy in redshift among streams makes them indistinguishable 
and hence compatible but inequivalent solutions result, as 
many as there are streams. The LAP method cannot dis- 
criminate among these solutions; multistreams indeed erase 
the memory of their past evolution. 

The second problem can only be overcome by casting aside 
part of the information contained in the sample and smooth- 
ing over the existing non-linearities to transform the multi- 
valued field into a single- valued one, typically with a smooth- 
ing length ~ 500 — 1000 kms -1 . The resulting smoothed 
field is clearly a less resolved representation of the underly- 
ing galaxy orbits, albeit the only tractable one. 

The advent of large galaxy redshift surveys strengthens 
the motivation to use the LAP method. Near all-sky red- 
shift surveys, e.g. PSCz, IRAS 1.2 Jy and ORS provide an 
excellent sky coverage (within a galactic latitude |6| <^ 8° for 
IR AS galaxies and |6] £ 20° for ORS), that may be extended 
further to cover the Zone of Avoidance via a Wiener recon- 
struction (Fisher et ai. 1995b; Zaroubi et al. 1999). They are 
therefore a fairly thorough representation of the underlying 
matter density field. Obviously the greater number of galax- 
ies in the sample the more accurate is the representation 
of the field, and this is best achieved with a redshift sur- 
vey. Real-space datasets require Tully-Fisher distance cal- 
ibrations of individual galaxies, and consequently the end 
result is a sparser sampling than is achieved with the same 
computational effort by measuring redshifts and angular co- 
ordinates. The goal of this paper is to exploit galaxy redshift 
surveys to the best effect and extract as much information 
from them as is possible; the main thesis put forward is the 
LAP method demonstrably breaks down the degeneracy in 
the determination of £7 m and b. This entails very tangible 
advantages. On the one hand, the freedom to investigate 
those two parameters separately permits us not to take the 
idea of bias seriously. A form of bias will certainly always be 
present in one form or another so that we can make sense 
of the galaxy number-count with respect to the underlying 
matter field. However, whether that is a linear or non-linear 
bias, the more one dissociates this phenomenological rela- 
tionship from our measurements of Q m , the more credible 
those measurements will be. This is indeed what LAP does. 
On the other hand, the LAP method produces a reconstruc- 



tion on the basis of the redshift-space sample alone, free of 
any proviso regarding the shape of the power spectrum. As- 
suming a given shape for P(k) unduly overconstrains the 
system, as will the addition of other datasets. 

In this article, I shall mainly apply the LAP method to 
the IRAS 1.2 Jy survey and study the predicted values of b 
and f2 m . The reconstructed IRAS 1.2 Jy velocity field is then 
compared with the Mark III velocity sample to seek a fine- 
tuning of the parameters. A more thorough undertaking, 
in terms of the quality of the sample, is to apply the LAP 
method to PSCz, which is by a factor of 3 a more densely 
sampled survey than IRAS 1.2 Jy, and it will be interesting 
to tackle this in future work. The article is structured as 
follows: Section 2 describes the LAP method in some detail 
and how to find solutions that are consistent with a redshift- 
space dataset; in Section 3 we test the method with several 
IRAS mock catalogues obtained via n-body simulations; in 
Section 4 we apply the method to the IRAS 1.2 Jy galaxy 
redshift survey, optimizing the predicted velocities with the 
Mark III dataset; finally, in Section 5 we summarize the 
main conclusions. 



2 THE LAP METHOD 

2.1 Redshift-space coordinates 

The redshift coordinates of galaxies are defined 



z = Hof+ f(f ■ v), 



(1) 



where r = (r, 8, ip) is the physical position, Ho is the present 
value of the Hubble parameter, v the peculiar velocity, and 
f a unit vector in the radial (line-of-sight) direction, z has 
units of velocity; its radial component is the redshift z r = cz, 
and the angular components are the same in both s-space 
and z-space, up to the distance scale. Henceforth we shall 
measure distances in km s , hence Ho is scaled out of the 
equations. In comoving coordinates, (ftl) reads 



S = X + x(x ■ Vjtt), 



(2) 



where the scale factor of the universe is normalized to 
a(to) = 1; ct(t, x) is the velocity potential, v = a -1 Va. Here- 
after we adopt to = 1. 



2.2 Dynamics 

The cosmological perturbations are derived from the action 



S = / dt dxC, 

J J sample 

where £ is given by 



3f2 m a 2 



(3) 



(4) 



5 is the density contrast and <j> the gravitational potential 
caused by the perturbations and 



C = 5+ - V- [(l + 8)v] 



(5) 



is the excess flux. The variations 8S/5vi = 5S /8cj> — yield 



v — -Va, 
a 



(6) 
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V 2 </> = — a 2 f2 m <5. 



(7) 



Similarly, SS/S5 — SS /5a = yield respectively 
£ = 0, (8) 
. IVa' 2 



2a 2 



+ = 0, 



(9) 



where we have eliminated v via (^) and we do not consider 
Oa. The field equations (|^),(^) are subject to the following 
boundary conditions: 

I Homogeneity of the density field at t — + 0. Density per- 
turbations grow from initial fluctuations of negligible ampli- 
tude: 



S(t ->(),£) « 0. 



(10) 



II Galaxy redshift survey at the present time. The galaxy 
number-count density p s in z-space constrains the real fields 
6~(x) and a(x) via 



p B (s) = x 



gals 



V 



1 + bS 
1 + a" 



(11) 



where the tilde denotes derivation along the radial direction, 
x is the radial comoving distance and b is the bias param- 
eter. Condition (7) is motivated by the CMB Sachs- Wolfe 
constraint 5 10~ 4 over r ~ 100, 000 km s _1 , so we accept 
that perturbations are negligible in the limit t —* 0. A proof 
for (II) is given in Appendix A. In order to solve (^|),(^|), we 
construct the trial fields: 



N 
N 

a = yjg n ft)a„(g) 

n=0 



(12) 



(13) 



where the basis functions f n ,g n are adjusted to numeri- 
cal convenience. SB94 considered f n = D(D — l) n , and 
g n = (D/D)f n , where D is the linear growth factor, nor- 
malized to unity at t = 1, so that the lowest-order series 
(p^),(p^|) are identical to the perturbative solutions. This 
is however strictly speaking not a compelling choice, and a 
sensible choice of orthogonal polynomials leads to an Ansatz 
of better convergence. As we have discussed in the Introduc- 
tion (point A), the sparseness of the dataset obscures the 
dynamical evolution and the LAP method is reduced to a 
numerical fit of the fields to the truncated equations, that we 
derive below, subject to ( |Io| ) , ([ll]) . In trying to approximate a 
function f(t) by orthogonal polynomials P m (t) in t 1, 
a weight function w(t) > tells us the relative importance of 
the errors spread over the domain. For a uniform w, f n are 
the [spherical] Legendre polynomials L m (t), whereas for a 
weight function that is larger at the endpoints ( |Io|) , (pd[) than 
throughout the trajectory, e.g. w(t) = (1 — t 2 )~ L/ ' ! (by shift- 
ing the domain from [0, 1] to [—1, 1]), the optimal choice are 
in this case Chebyshev polynomials T n (t). This choice mini- 
mizes the errors around the endpoints and it gives a greater 
weight to the solutions (matching the boundary conditions) 
in this region. In the analysis that follows, we shall adopt 
fn = T n and g n = a 2 f n . The fields 6„,ct n are expanded in 
terms of spherical harmonics, 



5n = £ Ji ( fcrX ) Y,m ' 

rim 

a n = ^2a< r " ) m ji(k r x)Yi m , 



(14) 
(15) 



where n are spherical Bessel functions. Substituting 

©,© into (|),(|) we get 



[ £a ' { Hmji(krx) + i(i A J} 



o ^ r(") 

--a 2 ft m ^ -^Tnji(k r x)Yi n 

rlmn 



(16) 



(17) 



the coefficients o/^™^ and Jj™^ are given in Appendix B. The 
boundary conditions (|l(]),(|ll]) then read 



ps = X 



N, 



gals 



V 



1 + bS(l,x) 



1 + a (l,x) 



(18) 



(19) 



where t is rescaled to the interval [—1,1] for convenience in 
using T n , and in @ we have used T n j-1) = (-1)™. The 
choice of basis functions of SB94 satisfy (jlq) by construction, 
and in our choice of basis functions the constraint is less 
trivial, but still it is easily tackled numerically. If we restrict 
ourselves to the interval < t < 1, then ( |l8| ) evaluated 
at t = eliminates all the Chebyshev polynomials of odd 
order. This is an equivalent approach but we shall adopt the 
convention above, —1 < t < 1. The constraint ([[9]) is the 
core of the problem as it is where all the information of the 
dataset is contained. The remainder of the paper will focus 
on the different ways one can use that constraint. 



2.3 Finding LAP solutions 

Substituting (p"2|)-(|l5"|) into equations (^),(j^), we get 



ji (k r x) Yi ri 



m 

N , 

) n — O rim ^ 



n— rim 
N 



(20) 



p,q— rim 
r' ' V ' m' 



+— 

x- 1 
and 

N 



xAJ}, q) ,(a) 

I'm 1 V / 



-fl k- 2 S {n) + ( — +2- \ a (n) 



T n ji(k r x)Yi 

■T, 



££[• 

n— rim 

N s 

= -\j2J2 T ' T *{ a 'rfrn<x'r%rn>MkrX)jl,(k r ,x) (21) 

i-i ri — fl rim. 



+ - 



p,q— rim 
T I'm' 



*A j£?(a) 



xA4S(a) 
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where the coefficients Ji^{8),J}^{a) are defined as in ( |B5[ ) 
in Appendix B and 8'^ m as in (|B3| ) via the trivial sub- 
stitution q — > 8. By multiplying (p0[),(|2l|) by T r jiYi m and 
integrating over all coordinates, we get 

JV JV 

^{T r T n )C S y 8 y n) +Y(T r T n )C y 'a v n) 



n=0 



n=0 



(22) 



P,<3=0 



n=0 



E 

n=0 



(T r (T„ 



•2-T n ))SJ4' 
a 



(n) 



JV 

^ ' (TrTpTq 
p,9=0 



9 1 V V ' 



(23) 



where y = (rim) and the angle brackets () for the Cheby- 
shev polynomials are defined in Appendix C. In deriving 
©,(§§), the coefficients C s y , C£, S s y , S%, D v ylyll and E y yly „ 
are calculated via Clebsch-Gordan coefficients for cross- 
products of Yim and via the standard orthogonality rela- 
tions for Yi m and ji, given in Appendix D. Cross-products 
of ji terms are estimated numerically. We proceed to solve 
(^),(^) numerically with the following iterative procedure. 
We first construct an Ansatz of the coefficients 6^ ,a4 that 
satisfies, to linear order, (p2;) , (p3|) as well as jl8|) , (Jl9|) - We 
start out with the galaxy number-count density p s . Follow- 
ing its definition in Appendix A, this quantity has units of 
inverse velocity, and we define its associated 2-space density 
contrast via 



4vriV, 



gals 



(1 + 8 S ), 



(24) 



where s ma x = cz max is the maximum redshift in the sample. 
Our first Ansatz entails b — 1 and linear evolution, so that 
6 S cx — V 2 a, and on inverting this relation to obtain the 
we estimate 5(x) tx 8 s (x + xa') by using the 



coefficients a y 



expression for the radial derivatives (B3). This yields a first 



Ansatz for 8 



derived from the dataset, that satisfies 



the linearized equations, given by the LHS of 



' C" c 4 " 




Oy 


s a s s 




8y 



(25) 



to do so, as a relationship of this kind between the redshift 
and real-space coordinates entails that we compare them via 
a Taylor expansion ji(k r s) « ji(k r x) + k r ct ' j[; an approxi- 
mation of this kind ~ O(0 2 ji) introduces an error of up to 
15% for I k, 10 as can be shown from (p32|) in Appendix B. 



2.4 Normal modes 

We have noted that the linearized equations ( p5| ) are a ho- 
mogeneous matrix system. If the determinant of the matrix 
is non-zero, then the only possible solution is S y — and 
d y = 0. We know however that ( p5| ) is also valid for linear 
fields, and these have non-vanishing coefficients. Therefore 
we conclude that the determinant of the system vanishes. 
Such a system of equations is tackled through the Singu- 
lar Value Decomposition (SVD) procedure. It factorises the 
singular matrix in (|2^) in a product of three matrices: two 
orthogonal matrices U and V, and a diagonal one W, which 
has one or more vanishing weights along the diagonal. After 
SVD, IH reads 



U 







Wl 



If 2 



V 



where the weights Wi, W2, 
Therefore, the vector 



Ny=V 



(26) 



. wn are non-zero real numbers. 



(27) 



gives a coordinate basis on which the first component, the 
normal mode, is unconstrained by the system (pHI). -/Vy is 
solely determined by (|l8|,(|is|). The rest of the components 
of N y (which are identically zero for linear fields) are func- 
tions of the normal mode. Therefore, one can rewrite the 
full non-linear system (p^),(p3|) in terms of the fields N y 
and this would be strictly speaking the natural basis to in- 
vestigate the underlying mode coupling induced by gravity. 
In the Fourier formulation with a set of basis functions like 
those used in SB94, /„ = D(D — l) n , it is easy to show 
numerically that the fc-th normal mode is given by 



l(0) 



<, 2 „,(°) 



(28) 



and (Sy) r = Sy 



(r) 



where the column vectors are (a y ); 
with r = 0, . . . , N. The solutions of the homogeneous system 
are then re-adjusted to satisfy ([l8]),([l9]) and we use these 
to construct the quadratic terms on the RHS of (p2j),(p3|). 
This leads to an inhomogeneous system that again we solve 
for 5 y ,d y . On each iteration we improve the solutions by 
least-squaring them to satisfy (]is|),(|l9|) to the best accuracy 
and we are also free to vary the parameters (6,f2 m ) for im- 
proved convergence. This procedure is very accurate, as we 
will show in the next sections, and it permits us to improve 
the estimate of the mapping x — > s at each iteration using 
the full non-linear relationship (^). At each iteration, the 
fields 8 y ,a y are used to obtain an estimate p s (s) of the RHS 
of (^). We then vary these fields to obtain a minimum of 
the quantity ^_(p s — p s ) 2 . Therefore we do not perform a 
jiYi m expansion of the dataset, and it is very convenient not 



This has a simple physical interpretation: (E8|) is a vanishing 
scalar for linear fields and thus its departure from zero gives 
us a measure of non-linearity. This quantity is determined 
by the boundary conditions. In the spherical harmonic for- 
mulation, the normal modes (equivalent to (Eg)) are 



iVf =Yh n {5^-k- 2 a^), 



n=0 



where 



dtw(t)D(t)T n 



(29) 



(30) 



where T] = 1 for n = and r\ = 2 otherwise. The quantity 
( [29] ) vanishes in the linear regime and, like (Ejj), its departure 
from zero is a measure of non-linearity. 
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2.5 Using the method in practice 



3 TEST OF THE METHOD 



The apparent mathematical complexity of the LAP method 
has precluded its wider use in practice. The fraction of pa- 
pers in the literature that employ LAP techniques to inves- 
tigate large-scale structure is minute in contrast to analyses 
based on perturbation theory techniques, such as POTENT, 
VELMOD and others. The latter unquestionably have the 
virtue of simplicity, and are as efficient as they are easy to 
implement. However, in practice the method described in 
this section entails no more complexity than programming 
an n-body code; an undertaking that merits the effort, so 
as to estimate b and f2 m , rather than merely j3. The chief 
difficulty resides in writing an algorithm for an effective nu- 
merical resolution of (|2^),(^). This may be a somewhat 
arduous task, but at any rate a very straightforward one 
with a very basic grasp of numerical methods. 

The LAP method is very flexible in its implementation. 
The basic input in the problem are the boundary conditions 
([i"o|), ( pd| ) and the procedure that is to be followed to find a 
stationary action linking both end-points is largely a matter 
of numerical convenience. The algorithm used in this section 
employs Chebyshev polynomials to fit the trial fields <5 and 
a to the dynamics. A myriad of other choices (e.g. binomial 
expansions, Legendre and Hermite polynomials, etc) is also 
feasible and thus the LAP implementation set out above is 
by no means a straight jacket recipe (for a more condensed 
presentation of the algorithm, see Susperregi 2000). 

In short, the algorithm can be summarized as follows. 

• A galaxy redshift survey is a dataset T> of points (z,<fi,9). 
Those raw data are transformed to a smoothed redshift- 
space field p s (s), given a smoothing length and a window 
function W(k). In this article we shall exclusively implement 
Gaussian smoothing. 

• The name of the game is to compute a fit for 8, a. The 
starting point is to make a linear Ansatz that is consis- 
tent with 5 S , which is derived from (^). This is achieved 
by inverting the relation 8 S cx — V 2 q and next estimating 
8 oc S s (x + xa'). 

• The linear Ansatz is the first input to be used in equa- 
tions (E3),(Ea). These yield the homogeneous system (p5|), 
which is our second port of call. The solutions 5 y ,a y obtained 
are least-square fitted to (|l^) , jli| ) . This requires adopting a 
value of b. 

• The adjusted values of S y ,a y are brought back to con- 
struct the RHS of (p2|),(p3|), and from there one computes 
the new S y ,a y in the LHS of (p2|) , (p3|) . This part of the oper- 
ation entails an assumed value for Q m . In the normal mode 
coordinates discussed in 2.4, the modes S y ,a y of the cosmic 
fields are merely excitation modes of a harmonic oscillator 
and the terms in the RHS of (p^),(p3|) represent nonlinear 
perturbations of those excitation modes. 

• Successive iterations of the procedure eventually yield 
the correct values of 5 y ,a y . The values of b and f2 m are read- 
justed in the process and their estimated values are those 
that result in the most rapid convergence of the solutions. 

The algorithm thus produces the cosmic fields and an 
estimate of the cosmological parameters. In the remainder 
of the article we shall investigate how to make the best use 
of the procedure and how to quantify the relative likelihood 
of different values of the cosmological parameters. 



We test the LAP method on mock catalogues derived from 
n-body simulations, using a Gaussian smoothing length of 
600 kms -1 . The IRAS 1.2 Jy power spectrum (Fisher et 
al. 1993) is adopted as a prior, and the simulations are 
performed over a periodic box L = 25,600 kms -1 with 
128 3 grid-points and 128 3 particles. The simulations are 
performed from Gaussian initial conditions, for the follow- 
ing values of the parameters: b = 0.8, 1.0, 1.2 and fi m = 
0.3,0.6,1.0. The fields are evolved forward in time until 
as « 0.7 over ~ 800 kms -1 , using a Gaussian cutoff. We 
choose a two-powerlaw functional form of selection function 
(Yahil et al. 1991): 



(r > r s 



ri + r' 



(31) 



r < r s ) — 1, where r s — 500 kms , r» = 5034 
, a = 0.483 and /3 = 1.79 (Fisher et al. (1995a); 



and 
km s 

we adopt the estimated central values of these parameters 
and will not test the fine detail of the variations of <j>(r) 
due to their errors) , and thus we compute the redshift-space 
dataset over a sphere of radius a; max ~ 17, 000 km s -1 . The 
resulting mock catalogue has an effective radius of ~ 13, 000 
km s -1 beyond which the galaxy number-count is sparse and 
is cut off for the purpose of the reconstruction. The num- 
ber of realizations are nine in total, and we denote d(b, fi m ) 
the 2-space mock samples derived in this way. Each dataset 
d(b, S7 m ) results from a unique pair of real-space fields 8, a 
which are the density contrast and velocity potential that 
we obtain via the n-body simulations. 

The tests are carried out by using d(b, f2 m ) as an input 
dataset in ( |l9| ) without any prior assumption on the real 
values of the parameters of the mock sample. We use ( jlil ) 
to solve (pj[),(p3|) following the iterative procedure given in 
§2.3 and derive the estimated fields 8, a for different values 
of the parameters 6,fi m . The likelihood of these parameters 
is estimated on the basis of the performance of the solutions 
8,a in terms of their convergence and ability to satisfy the 
constraints. We use a likelihood function given by the inverse 
of the x-squared sum of the differences of the fields between 
successive iterations in solving (E3),(p3|), i.e. 



s„ — s n 

as 



+ 



' Ctn-l 



(32) 



where n > 25, as ~ 0.20, a a is a normalization factor that 
rescales the coefficients a n so that a becomes a dimension- 
less field within the range — 1 Si a Si 1 and we have used 
N = 10 and / < 15 and an initial linear Ansatz. The results 
are shown in Fig. 2. The likelihood contours are the LAP 
reconstructions and the crosses on all nine panels of Fig. 2 in- 
dicate the values of the real parameters in each mock dataset 
on the (6, fi m ) plane. As can be observed, the likelihood con- 
tours are certainly well correlated with the location of the 
crosses, where the innermost contours mark the level of 95% 
likelihood, that in all cases lie in the neighbourhood of the 
real values of the parameters. The likelihood contours show 
an approximately elliptical shape, with the major semiaxes 
tilted at approximately 45 degrees, suggesting a correlation 
between both parameters that merely arises in the numerical 
computation. The estimates in the reconstruction are fairly 
good, with a trend in underestimating slightly the values of 
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Figure 2. Likelihood contours for the reconstruction of the nine datasets d(b, C m ). The cross on each panel indicates the real values of 
(b, flm) in each reconstruction, and the likelihood contours are computed following (B2J) with a suitable normalization. The concentric 
contours represent a likelihood of 95%, 75%, 50%, 25% and 10% from the inner curves to the outer, on the two upper rows, and 95%, 
75% and 50% on the lower row. 



both parameters. The best reconstructions are for the inter- 
mediate value of the density parameter fl m = 0.6, shown in 
the second row. In this case the crosses actually lie within the 
highest likelihood contours, with very little scatter. Overall, 
in the nine reconstructions the rms scatter in b and Q m lie 
within the region 0.26 £ cr£ £ 0.44, 0.15 ?S a\ ^ 0.32. The 
largest scatter in f2 m is for f2 m = 0.6, and a similar situation 
arises with b, whereby the intermediate value b — 1.0 has the 
larger error. 

The effect of underestimating the true values of the pa- 
rameters is systematic and can be calibrated. This effect can 
be largely ascribed to the unconventional choice of likelihood 
estimator (|32|). One could argue that, for slowly varying vari- 
ances, A oc b~ 2 (chiefly from the S part of the RHS of @) 
and therefore smaller values of the bias factor (and conse- 
quently, by correlation, also of fl m ) are favoured. However, 
it is not straightforward to disentangle the dependence of 
the solutions on the parameters after successive iterations. 
The likelihood estimator used is thus to some extent biased. 
However, we find that the criterion of convergence given by 
the RHS of (^j|), suitably normalised, is the sharpest dis- 
criminator to pin down the best estimates of the cosmolog- 
ical parameters. We have carried out numerous tests with 



more conventional likelihood estimators (e.g. Fisher likeli- 
hood matrix, etc) obtaining much poorer results than with 

©• 

Fig. 3 shows the density constrast reconstructions for 
the same datasets d(b, fi m ). The reconstructed density con- 
trast 5lap is shown on the horizontal axis plotted point- 
by-point within the selected spherical volume (r ~ 13, 000 
kms -1 ) against the real density contrast of the mock 
datasets. A solid line of slope 1.0 is plotted across each panel 
that does not correspond to the regression line on each panel 
though the differences are tiny. The slopes of the regression 
lines lie within the range 0.99 ± 0.08. The rms value cor- 
responding to the random and numerical errors lies in the 
range 0.19 <> as ?S 0.28. The reconstructions in Fig. 3 have 
been carried out with a prior knowledge of the values of b, fi m 
for each dataset. Alternatively, the test can be carried out by 
putting together the procedure followed to obtain the like- 
lihood in Fig. 2 and investigate the scatter resulting in the 
plots 6 m ock vs - Slap for different values of 6,f2 m . Supposedly 
estimating the values of 6, fi m and finding the optimal cor- 
relation between 5 m0 ck,<5LAP ought to be two not unrelated 
operations. However these two appear to be fairly indepen- 
dent: it turns out that whereas (^) gives us the correct 
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likelihood estimates following the criterion of convergence 
of the solutions at each iteration, the variations in as for a 
large range of b,Sl m are fairly small, and as (as computed 
from tests such as the nine reconstructions in Fig. 3) is too 
insensitive to be helpful in the estimate of the parameters. 
Therefore the tests show that the estimate of the parameters 
and the reconstruction of the fields are two operations that 
are to a large extent independent. For an arbitrary sample, 
one would thus first compute (^), pick the values of b, Q m 
at the maximum of the likelihood surface and use these to 
solve the equations to compute 8, a. Similarly, Fig. 4 shows 
the comparison of the LAP results with the mock data in 
the reconstruction of the velocity potential. The values of 
the fields have been scaled to a max and are dimensionless. 
It is apparent that the regression line is in all cases slightly 
greater than unity, with a more accentuated tilt for larger 
values of (fo,S7 m ). The smaller values of a adjust better to a 
slope of unity, but with larger scatter than larger a. 

Fig. 5 shows a cross-section on the Z = plane of a 
particular velocity field reconstruction, that of the dataset 
d(b — 1.0, Sl m = 0.3). The figure shows several prominent 
features of the underlying density field in this case: three 
overdense regions to which the field vectors converge, on the 
lower left, middle right and upper left parts of the panel, 
and two prominent underdense regions, from which the ve- 
locities diverge, one at the central region and another one at 
the middle-left boundary of the circle. It is apparent that the 
LAP velocities are not vanishing in the normal direction of 
the boundary surface of the selected subvolume, and there- 
fore the customary Neumann spatial boundary conditions 
employed on spherical Bessel functions (i.e. vanishing nor- 
mal velocities at the boundary) do not apply. We note that 
spatial boundary conditions are unnecessary in the LAP re- 
construction, thus we have not brought up the issue in §2. 
The velocity field agrees within 10% accuracy with the n- 
body exact field within 78% of the selected volume, and 
the remaining 22% differs from the mock sample velocities 
by an error of <* 10% (shown in Fig. 4 by the regions en- 
closed by the solid curves) and withing this volume 6% dif- 
fers by an error <^ 20% (regions enclosed by broken curves). 
These regions are mostly located in the neighbourhood of 
peaks, right at the very slopes, where the largest velocities 
are found. The central regions of peaks and troughs are very 
accurately reconstructed, and it is indeed the intermediate 
regions that yield 8 points with greater scatter in Fig. 3 and 
worse velocity reconstructions in Fig. 5. 



4 BIAS AND Q m FROM IRAS 1.2 Jy 

We apply the LAP method to the IRAS 1.2 Jy sample 
(Strauss et al. 1990,1992; Fisher et al. 1995a) in the same 
way as we have used it in the reconstruction of the mock 
catalogs in §3. IRAS 1.2 Jy is not the largest existing near 
all-sky galaxy redshift catalogue, and it is now superseded 
by PSCz (Canavezes et al. 1998) which contains ~ 15, 000 
galaxies, so this application is simply an illustration on how 
the LAP method can be used to break the degeneracy in the 
estimates of b and f2 m . Other large redshift samples of par- 
tial coverage can also be looked at with the LAP method, 
e.g., Las Campanas and the forthcoming Anglo- Australian 
2dF (~ 250, 000 galaxies) and US Sloan Digital Sky Survey 



12 




X 

(10 3 kms" 1 ) 

Figure 5. Z = plane reconstruction of the velocity field within 
a selected subvolume x max ^ 12, 000 km s — 1 for the mock sample 
d(b = 1.0, Q m = 0.3). The solid lines enclose regions where the 
reconstruction entails an error |# moc k — «l AP I /"mock ^ 0.10 and 
within the broken lines this error is 0.20. 



(SDSS) (~ 10 6 galaxies and 25% coverage), with the caveat 
that boundary regions will be a source of propagating er- 
rors in the dynamical evolution. Even so, a large number 
of galaxies in a redshift survey of limited coverage can pro- 
vide a good representation of the underlying density field, 
almost definitely outweighing the disadvantages of sampling 
a partial region of the sky, and it will be thus predictably 
worthwhile to apply the LAP method to those surveys. The 
IRAS 1.2 Jy sample contains 5320 galaxies distributed over 
87.6% of the projected celestial sphere. The remaining un- 
sampled 2.4% is an approximately disk-shaped region at a 
galactic latitude |6| 5s 5°. 

We adopt a Gaussian smoothing length of 1200 km s _1 , 
and make no assumption regarding the power spectrum. The 
data diRAS are distributed within a spherical region of radius 
i ma j ~ 15, 000 km s^ 1 . We use the dataset in a similar fash- 
ion as the mock samples d(b, f2 m ) in the previous section to 
derive the a;-space fields 8, a. In §3 we have established that 
as, a a are fairly insensitive to the values of b, f2 m . One can 
thus set out to investigate the likelihood function X(b, S7 m ) as 
defined in ( pS2] ) prior to determining the reconstructed fields. 
Evidently this is the simplest way to proceed for, unlike in 
§3, we do not have any clue about the real-space underlying 
fields (such as <5 moc k, a m ock in §3) to compare them with the 
reconstructed fields. 

The likelihood contour plot is shown in Fig. 6. Clearly 
the largest values of the likelihood function are centered 
around b ~ 1 and small Sl m . From the test of the LAP 
method in §3 with n-body simulations we already know that 
the likelihood function (K2|) underestimates both b and fi m , 
as is apparent in all nine panels of Fig. 2. We accept this 
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Figure 3. Density field reconstructions of the nine datasets d(b, f2 m ). The smoothed density contrast of the mock samples (vertical axis) 
is compared at each point within a selected spherical volume of the 128 3 grid with the LAP-reconstructed densities (horizontal axis) over 
a sphere of radius ~ 13, 000 km s . The systematic errors are caused by the sparseness of the sampling. 



trend is fairly inherent to the numerical application of the 
method and thus infer that the result presented in Fig. 6 
is no different in this respect, and therefore the real values 
of the parameters are situated somewhat above their max- 
ima in the likelihood function. From Fig. 2 one can quantify 
these errors to be of the order of Afi m ~ 0.12, Ab » 0.15. 
Therefore, we infer that in Fig. 6 the likelihood maxima and 
the real values of the parameters are likely to be offset by a 
similar margin of error. At face value, Fig. 6 estimates that 
the most likely values of the parameters are fl m « 0.18 and 
b ~ 0.94. If we offset these estimates by the errors derived 
from Fig. 2, then the likely "real" values of the parameters 
that we obtain for Fig. 6 are fi m « 0.31 and b » 1.1. As a 
matter of fact, these offset values are still within the region 
enclosed by the 95% confidence contour. 

To put our results in perspective with previous analyses 
of IRAS 1.2 Jy, we have overlaid on the contour plot of Fig. 6 
two previous estimates of (3 = Q^ 6 jb. An estimate by Willick 
et al. (1997a) yields fa = 0.49 ±0.07 (shaded region A) and 
an estimate by Sigad et al. (1998) yields fa = 0.89 ± 0.12 
(shaded region B). The estimate of Willick et al. (1997a) is 
clearly in better agreement with our results as the location 
of the offset maximum of the likelihood is contained within 
the shaded region A that corresponds to the error margin 



of their estimate. The estimate given by shaded region B 
is consistent with a scenario b ps 1, f2 m « 1, which in our 
analysis falls well outside the 10% likelihood contour. 

Fig. 7 shows a z-space comparison between the recon- 
structed fields and the dataset. The data on the horizon- 
tal axis, <5£api i s obtained from the reconstructed x-space 
fields 8, a via (pi]). The combination of both fields via the 
relationship 5lap(i) oc S a (x + xoe' LAP ) permits us to recon- 
struct 8 a which is our only possible point of comparison with 
Siras, and this is shown in Fig. 7. The vertical axis shows 
the 2-space data points of the smoothed IRAS 1.2 Jy sam- 
ple. The data are plotted in a point-by-point comparison for 
all the grid points within the selected subvolume. A solid 
line of slope 1.0 is plotted across the diagonal of the plot. 
The slope of the regression line is slightly over the diagonal 
line, at approximately 1.03. The corresponding rms due to 
random and numerical errors in the LAP reconstruction is 
o ~ 0.27. The values of the parameters that have been used 
in the reconstruction are b — 1.0, fi m = 0.3. 

4.1 Velocity fields 

The resulting velocity field for the parameters of Fig. 7 
is shown in Fig. 8. The six panels show the reconstructed 
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Figure 4. Velocity potential field reconstructions of nine datasets d(b, C m ). The velocity potential values are scaled to o ma x; so that 
they are dimensionless and consigned to the range —1.0 ^ a ^ 1.0. The smoothed velocity potential of the mock samples (vertical axis) 
is compared at each point within the same selected volume as in Fig. 3. 



IRAS 1.2 Jy fields Slap and vlap in supergalactic coordi- 
nates, for three different slices Z — 0, ±2000 km s . The 
velocity panels on the right column correspond to the den- 
sities on the left, at the same value of Z. The velocity field 
follows the main features observed on the <5lap field, with a 
general flow towards the overdense regions and outflow from 
voids. The largest velocities are located in the intervening 
regions between overdense and underdense regions, e.g. in 
Z = (middle panels), large infall velocities are visible in 
the vicinity of the Comma supercluster (0,80,0), the Hydra- 
Centaurus (H-C) supercluster (-30,15,0), and Perseus-Pisces 
(P-P) (50,-5,0). In Z — the largest velocities are located 
at the lower right region of the H-C overdensity maximum, 
and also to the left of the P-P maximum. There is a veloc- 
ity flow from the main void on the lower left of the figure, 
in the direction of Virgo, and it splits up to left and right, 
in manner of a ridge, to create an outflow in opposite di- 
rections, towards H-C and P-P. In the case of Z = —2000 
km s -1 (lower row), large velocities are also present around 
the steeper regions of the prominent overdensities, follow- 
ing a similar pattern as in Z — 0, whereas the field shows 
more erratic features in Z — 2000 km s _1 (upper row), where 
the outflow from the main void (centre left) shows a general 



trend towards the main overdense features but is at the same 
time prone to local variations. 

The results presented in Figs. 6-8, can be optimized 
by using the Mark III velocity redshift survey to pin down 
b,flm more accurately. We shall pursue this and look for the 
optimal values of b,£l m by computing the LAP solutions that 
satisfy 

S ^(«LAP — "Marklll) 2 = 0, (33) 

where S denotes a variation, not the density contrast. In 
practice, this is achieved as follows. One adds (^) to the two 
already existing constraints of the LAP method (|l8|),(|l9[). 
Those are tackled in the manner summarized in §2.5. In 
actual terms, it's far more practical to deal with (^) in 
terms of the velocity potential, so what we have done in the 
present analysis is in reality to compute QMarH n fr om the 
smoothed observed velocity field, and thus used (p3|) in the 
manner of a second constraint on a. 

The comparison with the WMarklll data sets further con- 
straints on the likelihood contours of Fig. 6 as is shown be- 
low. Mark III contains approximately 3,400 galaxies, which 
are compiled from several sets of elliptical and SO galax- 
ies (Willick et ai. 1995,1996,1997a). The sample spans out 
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Figure 8. <5lap and i?LAP fields for IRAS 1.2 Jy. SGX and SGY units are in 100 kms -1 , spanning over a sphere of radius ~ 8000 
km s~ 1 . Left: from top to bottom panels, density contrast for a Gaussian smoothing of 1200 km s — 1 , for Z = 2000, 0, —2000 km s _1 . 
Thick solid line corresponds to <5 = 0, continuous contours are <5 > and slashed contours are <5 < 0; contour spacing is 0.2. Right: from 
top to bottom, reconstructed velocities at same values of Z. 



to ~ 6000 kms -1 , though in some directions it is irreg- 
ularly sampled to ;r m ax ~ 8000 km s _1 and Xmin ~ 4000 
km s _1 . The distances are inferred via forward Tully-Fisher 
and D n — a distance indicators which may entail an er- 
ror in the region 17-21%. Mark III predicts a bulk flow 



vb ~ 194 ± 32 km s _1 towards the Shapley concentration 
(Zaroubi, Hoffman & Dekel 1999) (for a low- resolution Gaus- 
sian smoothing ~ 1200 km s - , within a sphere r ~ 6000 
km s _1 ), in contrast to vb ~ 250 — 400 km s _1 that is esti- 
mated in most other samples, including PSCz (a compilation 
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Figure 6. Likelihood contours for IRAS 1.2 Jy. The concentric 
contours represent a likelihood of 95%, 75 %, 50%, 25 %, 10% from 
the inner curve to the outer. The shaded region A corresponds to 
the estimate of by Willick et al. (1997a) and the shaded region 
B corresponds to an estimate of /3 by Sigad. et al. (1998). 



Figure 9. Solid contours represent the likelihood for IRAS 1.2 
Jy as in Fig. 6, and dotted contours represent the likelihood in 
the IRAS 1.2 Jy/Mark III comparison following (^)- The relative 
likelihood of the concentric contours is as in Fig. 6 in both solid 
and dotted. 
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Figure 7. Redshift-space density contrast in the LAP recon- 
struction versus the corresponding IRAS 1.2 Jy data for a Gaus- 
sian smoothing of 1200 km s -1 within a spherical region of radius 
12,000 km s -1 . The field <5 LA p is evaluated at 6 = 1.0, 
tt m = 0.3. 

of vb estimates is summarized in Dekel 1999b). Siras and 
^Markin are consistent with mildly non-linear gravitational 
instability and linear bias (Sigad et al. 1998), though there 
are some differences, e.g. the Mark III sample appears to 
show a strong shear across the Hydra-Centaurus complex 
that is absent in IRAS 1.2 Jy (as indeed also in ORS). Re- 
cent papers have studied in detail the differences between 
the IRAS 1.2 Jy and Mark III velocity and density fields 
(Sigad et al. 1998; also Dekel et al. 1999 following an im- 
proved version of POTENT). 

We consider the Mark III sample with a Gaussian 
smoothing length of 1200 kms -1 . The data are carefully 



corrected for Malmquist biases (following the recipe set out 
in Sigad et al. (1998) for the preparation of the data), and 
the distances of 1,241 objects are modified as a result. The 
LAP method is solved for IRAS 1.2 Jy within spherical vol- 
ume of radius a; max ~ 15, 000 km s" 1 , and the minimization 
fit with Mark III ( |33| ) is done within a spherical subvolume 
of radius (x) ~ 6000 km s _1 . Therefore most of the volume 
of the LAP solutions remains free of the constraint ( [33| ) and 
the fraction of the volume where «lap is least-squared to 
^Markin is only 0.064. Naturally such a small fraction fore- 
casts an almost negligible impact in the fine-tuning of the 
parameters, unless the fields differred drastically to start 
with, which they do not. The «lap solution in the remain- 
der of the volume is indirectly affected by this fit, and the 
variations in modulus Avlap outside the comparison sub- 
volume are Ss 12%. 

Fig. 9 shows the likelihood contours for (b,Q m ) com- 
puted via the adjustment entailed in (^). The solid con- 
tours are the purely IRAS 1.2 Jy prediction, as in Fig. 6, 
and the dotted contours are the result of the comparison 
with Mark III. The contours are ever so slightly shifted to- 
wards greater values of the parameters and, as expected, the 
effect is small. The shift towards larger b,S} m is not in fact 
an altogether undesirable modification, as we have already 
discussed that the LAP solutions are found to be per se off- 
set to smaller values than their "real" values. The important 
conclusion to be drawn from Fig. 9 is that the comparison 
with Mark III is entirely consistent with the predictions for 
b and f2 m extracted from the IRAS 1.2 Jy sample alone. 



5 CONCLUSIONS 

The LAP method provides a practical means to break the 
degeneracy between Q m and b in galaxy redshift surveys. 
The method is employed in the manner of a nonlinear con- 
straint on the redshift-space dataset and, although in for- 
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mulation it comes across as algebraically cumbersome, it is 
of considerable simplicity and efficiency from the numerical 
point of view. The method is sound in that it does not re- 
quire an a priori approximation of the map x — > s to pin 
down the solution and it provides considerable freedom to 
ascribe relative importance to the data available, i.e. the 
initial and final endpoints, to which we wish to invariably 
assign greater weight than intermedate stages of which little 
or no data are available. 

The method can prove significant to measure fi m in 
the latest largest samples, and extract the most accurate 
information prior to comparison with other datasets, such 
as the CMB radiation power spectrum and SN data. One 
important challenge for the future is to attain a better grasp 
of the concept of bias and this will be probably achieved via 
microlensing data and n-body simulations of the formation 
of galaxies and clusters from primordial fluctuations, rather 
than from galaxy redshift surveys. Once a model of bias is 
adopted on a sound footing, then clearly the LAP model is 
impeccable in producing an estimate of f2 m . In the simple 
linear bias model we have employed we have totally relegated 
any consideration of scale-dependence in b. This is a point I 
have deliberately omitted for simplicity. Thus, the estimates 
computed in this paper ought to be regarded qualitatively 
as weighted averages of the "real" b over different scales, if 
indeed scale-dependent bias models are to be believed. 

In this paper, we have employed the likelihood func- 
tion (^2|) to investigate the values of b,Q. m . Clearly this is 
not a unique choice. However, our choice is guided by the 
argument of relative convergence of the solutions, which is 
justifiably a reasonable criterion to get close to the "real" 
solutions. In view of the performance of the A function in 
the reconstruction of the mock samples, this choice does not 
appear to be totally off the mark. A potential reason for con- 
cern could be the offset observed between the maxima of the 
likelihood functions and the real values of the parameters in 
the n-body simulations. However the recurrence of this off- 
set in a predictable manner lends strength to the argument 
that it arises as a numerical fault that is easy to account 
for systematically in the analysis of the datasets. The re- 
constructions of the fields are, on the other hand, of consid- 
erable accuracy and no numerical defficiency or hindrance 
is observed. The application of the method to IRAS 1.2 
Jy predicts the parameters to be fairly accurately located 
in the immediate neighbourhood of the maxima Q m ~ 0.3 
and b m l.f, which is found to be most compatible with the 
estimate of j3 given by Willick et al. (1997a). In a flat uni- 
verse such predicted values are perfectly consistent with a 
non- vanishing cosmological constant or a quintessence scalar 
field component. The likelihood examined in this way is only 
very slightly modified when the velocities predicted via the 
LAP method are finely-tuned with data from the Mark III 
sample. The shift of the predicted values is towards slightly 
greater values of the parameters but it remains comfortably 
consistent with the results obtained from IRAS 1.2 Jy. 
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APPENDIX A: ORBIT-CROSSING IN 
REDSHIFT SPACE 

We shall prove the boundary condition (0). The number- 
counts of galaxies n in x-space and z-space satisfy, by con- 
servation of the number of galaxies: 



dn(s) = dn(xi), 



(Al) 



for all streams at the same redshift, s = Xi + x(x ■ V^a,). 
In our analysis we shall only consider single-valued solu- 
tions, and therefore there is just one stream only in (Al), 
i.e. dn(s) = dn(x). Hence 



p s (s) dQ 



dn(s) 
ds 



. b5) N^sdx 
V ds 



(A2) 



where n(s) is the galaxy number-count, dtt a solid angle 
element and the x-space selected volume of the sample is 
V ~ fTTX^ax, and 

s = x ■ s = x + a' . 

Therefore 
da; _ 1 
ds - 1 + a" ' 



(A3) 



(A4) 



and substituting this in ( A.2 ) , we get 



iVgals ( 1 + bS 



V 



l + a' 



(A5) 



In the case of multistreams, the RHS of (A5) is integrated 



over all streams, bearing in mind that turn-around regions, 
which occur at 8 3> 1 and for which ds/dx = 0, are ex- 
cluded from the sum. An example of such a region in shown 
in Fig. Al. An initial saddle-point ds/dx = on the s(x) 
curve starts the creation of a turn-around region. At the 
stage shown in Fig. Al, both points A and B satisfy this con- 
dition and obviously they departed from an initial saddle- 
point A = B. The region spanning between A and B is 
three- valued (each redshift in the interval zb < z < za 
corresponds to three x positions), whereas za and zb are 
bivalued. To make such an scenario tractable, we need to 
replace s(x) over the interval zb < z < za by a monotonic 




A 



B 



x 



Figure Al. Illustration of a turn-around region. 



curve that matches the existing curve at zb and za and its 
first derivative. This is obviously tantamount to applying a 
larger smoothing length than the existing one to erase the 
overdense region that is the cause of the turn-around. 



APPENDIX B: EVALUATION OF RADIAL 
DERIVATIVES 

The radial derivative of the velocity potential coefficients 
( p^| ) can be written as 

d 



(i ~n = ^2a' < r "l l ji(k r x)Yi 

( 

ji(u) = (21 + I)" 1 Iji-^u) -(1 + l)i,+i(u) 



where, using the equality 

_d 

du 

we have 



rv' (n) - k 



Similarly 



I 



+ 1) (n) 

(2Z + 3) (21 - i)"r(i-i) m 



(Bl) 



(B2) 



(B3) 



//(«) ,2 f (1 + 2) (n) 

" S (21 + 3) (2Z + 5) r ('+ 2 ) m 



(« + l) 2 



r 



(2l + 3)(2l + l) (2Z-1)(2Z + 1) 



„(«) 



(B4) 



(2Z - 1) (2Z - 3) 
On the other hand, the coefficients J}^ given in ( [l6|) are 



ft") 



E 



a(l,m + 1) {n) 



W(m+1) ( lXl ~ X V 
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(3(1, m — 1) („) . . („) „ 

O rl(m— 1) V ' rim 



where 



a(£, m) 



/(/ + 1) - m(m - 1) 



(3(1, m)= 1(1 + 1) -m(m+l) 



1/2 



1/2 



(B5) 

(B6) 
(B7) 



APPENDIX C: CHEBYSHEV POLYNOMIALS 

The Chebyshev polynomials are defined T^cos^) = cos(n#) 
(following the normalization of Abramowitz & Stegun 1972) . 
We define the angle brackets {, ) according to the orthogo- 
nality properties of T n (e.g. Courant & Hilbert 1989): 



(u) 



dtw(t)u(t), 



(CI) 



where w(t) = (1 — t 2 ) is a weight function and therefore 



(C2) 



for n/0 and (Tq ) = tt. In (|22)(^) we encounter two types 
of angle brackets to evaluate (other than (|C2j) : (T n T m ) and 
(T n T, n T r ) (we have deliberately omitted (n m T n T m ), by ap- 
proximating fi m by a constant, and ditto for H . The second 
type of product is trivially transformed into (|C2|) via 



(C3) 



for n > m, and the first requires a little numerical manipu- 
lation using the relation 



(1 - r)f n 



-ntT n + T n -i. 



(C4) 



APPENDIX D: ORTHOGONALITY 
RELATIONS 

The orthogonality relations for the spherical harmonics and 
the Bessel functions are respectively 



dap / d(cos 0)Y im y; 

l 



dxx ji(k r x)ji(k 3 x) 



j^k^+xj^krx) 



(Dl) 



<5 rs .(D2) 
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